1 Wymagania

1.2 Executive summary

raport powinien zaczynać się od rozdziału podsumowującego całą analizę, streszczającego najważniejsze spostrzeżenia analityka

1.3 Lista wymagan minimalnych [kolejnosc dowolna]

Kod wyliczający wykorzystane biblioteki. Kod zapewniający powtarzalność wyników przy każdym uruchomieniu raportu na tych samych danych. Kod pozwalający wczytać dane z pliku. Kod czyszczący dane (np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych). Sekcję podsumowującą rozmiar zbioru i podstawowe statystyki. Szczegółową analizę wartości atrybutów. Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji. Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie. Sekcję próbującą stworzyć klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji. Analizę ważności atrybutów najlepszego znalezionego modelu. Dodatkowo punktowane będzie wykonanie analizy typowej dla danych klinicznych, np. regresji logistycznej wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors) lub regresji Coxa (ang. Cox Proportional-Hazards Model).

1.4 Dodatkowe uwagi

Analityk nie musi, a nawet nie powinien, ograniczać się do powyższych punktów. Wszelkie dodatkowe techniki analizy danych, wizualizacje, spostrzeżenia będą pozytywnie wpływały na ocenę.

Ewentualne konkluzje, znalezione zależności warto potwierdzić dokonując sprawdzenia istniejących wyników badań w literaturze naukowej (np. na Google Scholar czy PubMed).

2 Wykonanie

2.1 Podsumowanie

TBA

2.2 Import bibliotek

library(xlsx)
library(DT)
library(knitr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(janitor)
## 
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(imputeTS)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lares)
## 
## Attaching package: 'lares'
## The following object is masked from 'package:janitor':
## 
##     crosstab
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

2.3 Wczytanie i wstepna analiza danych

2.3.1 Wczytanie danych,

filename <- 'res/wuhan_blood_sample_data_Jan_Feb_2020.xlsx'
raw_data <- read.xlsx(filename, 1)
raw_data <- as_tibble(raw_data)

2.3.2 Sprawdzenie surowych danych

# OK
#how many NA in variables
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770596
# var min i max
min(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)
## [1] -1
max(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)
## [1] 50000
glimpse(raw_data)
## Rows: 6,120
## Columns: 81
## $ PATIENT_ID                                                    <dbl> 1, NA...
## $ RE_DATE                                                       <dttm> 2020...
## $ age                                                           <dbl> 73, 7...
## $ gender                                                        <dbl> 1, 1,...
## $ Admission.time                                                <dttm> 2020...
## $ Discharge.time                                                <dttm> 2020...
## $ outcome                                                       <dbl> 0, 0,...
## $ Hypersensitive.cardiac.troponinI                              <dbl> NA, N...
## $ hemoglobin                                                    <dbl> NA, 1...
## $ Serum.chloride                                                <dbl> NA, N...
## $ Prothrombin.time                                              <dbl> NA, N...
## $ procalcitonin                                                 <dbl> NA, N...
## $ eosinophils...                                                <dbl> NA, 0...
## $ Interleukin.2.receptor                                        <dbl> NA, N...
## $ Alkaline.phosphatase                                          <dbl> NA, N...
## $ albumin                                                       <dbl> NA, N...
## $ basophil...                                                   <dbl> NA, 0...
## $ Interleukin.10                                                <dbl> NA, N...
## $ Total.bilirubin                                               <dbl> NA, N...
## $ Platelet.count                                                <dbl> NA, 1...
## $ monocytes...                                                  <dbl> NA, 1...
## $ antithrombin                                                  <dbl> NA, N...
## $ Interleukin.8                                                 <dbl> NA, N...
## $ indirect.bilirubin                                            <dbl> NA, N...
## $ Red.blood.cell.distribution.width.                            <dbl> NA, 1...
## $ neutrophils...                                                <dbl> NA, 6...
## $ total.protein                                                 <dbl> NA, N...
## $ Quantification.of.Treponema.pallidum.antibodies               <dbl> NA, N...
## $ Prothrombin.activity                                          <dbl> NA, N...
## $ HBsAg                                                         <dbl> NA, N...
## $ mean.corpuscular.volume                                       <dbl> NA, 9...
## $ hematocrit                                                    <dbl> NA, 3...
## $ White.blood.cell.count                                        <dbl> NA, 3...
## $ Tumor.necrosis.factorÎ.                                       <dbl> NA, N...
## $ mean.corpuscular.hemoglobin.concentration                     <dbl> NA, 3...
## $ fibrinogen                                                    <dbl> NA, N...
## $ Interleukin.1β                                               <dbl> NA, N...
## $ Urea                                                          <dbl> NA, N...
## $ lymphocyte.count                                              <dbl> NA, 0...
## $ PH.value                                                      <dbl> 7.415...
## $ Red.blood.cell.count                                          <dbl> NA, 4...
## $ Eosinophil.count                                              <dbl> NA, 0...
## $ Corrected.calcium                                             <dbl> NA, N...
## $ Serum.potassium                                               <dbl> NA, N...
## $ glucose                                                       <dbl> NA, N...
## $ neutrophils.count                                             <dbl> NA, 2...
## $ Direct.bilirubin                                              <dbl> NA, N...
## $ Mean.platelet.volume                                          <dbl> NA, 1...
## $ ferritin                                                      <dbl> NA, N...
## $ RBC.distribution.width.SD                                     <dbl> NA, 4...
## $ Thrombin.time                                                 <dbl> NA, N...
## $ X...lymphocyte                                                <dbl> NA, 2...
## $ HCV.antibody.quantification                                   <dbl> NA, N...
## $ D.D.dimer                                                     <dbl> NA, N...
## $ Total.cholesterol                                             <dbl> NA, N...
## $ aspartate.aminotransferase                                    <dbl> NA, N...
## $ Uric.acid                                                     <dbl> NA, N...
## $ HCO3.                                                         <dbl> NA, N...
## $ calcium                                                       <dbl> NA, N...
## $ Amino.terminal.brain.natriuretic.peptide.precursor.NT.proBNP. <dbl> NA, N...
## $ Lactate.dehydrogenase                                         <dbl> NA, N...
## $ platelet.large.cell.ratio.                                    <dbl> NA, 3...
## $ Interleukin.6                                                 <dbl> NA, N...
## $ Fibrin.degradation.products                                   <dbl> NA, N...
## $ monocytes.count                                               <dbl> NA, 0...
## $ PLT.distribution.width                                        <dbl> NA, 1...
## $ globulin                                                      <dbl> NA, N...
## $ γ.glutamyl.transpeptidase                                    <dbl> NA, N...
## $ International.standard.ratio                                  <dbl> NA, N...
## $ basophil.count...                                             <dbl> NA, 0...
## $ X2019.nCoV.nucleic.acid.detection                             <dbl> NA, N...
## $ mean.corpuscular.hemoglobin.                                  <dbl> NA, 3...
## $ Activation.of.partial.thromboplastin.time                     <dbl> NA, N...
## $ High.sensitivity.C.reactive.protein                           <dbl> NA, N...
## $ HIV.antibody.quantification                                   <dbl> NA, N...
## $ serum.sodium                                                  <dbl> NA, N...
## $ thrombocytocrit                                               <dbl> NA, 0...
## $ ESR                                                           <dbl> NA, N...
## $ glutamic.pyruvic.transaminase                                 <dbl> NA, N...
## $ eGFR                                                          <dbl> NA, N...
## $ creatinine                                                    <dbl> NA, N...
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770596

2.4 Transformacja danych

(np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych)

# NOT OK
#ogarnac te brakujace daty

#replace -1 with NA
raw_data[raw_data==-1]<-NA

#PATIENT_ID
id_filled <- raw_data %>% fill(PATIENT_ID)


#OTHER NA VALUES
#is there any rows that doesn't have NA value?
full_row_count <-dim(id_filled[complete.cases(id_filled),])[1]
print(full_row_count)
## [1] 0
# NA values in variables:
mean(is.na(id_filled %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770748
# there are two patients with negative days, as their only blood sample results arrived one day after their clinical outcome


#remove rows where all variables are empty
vars <- colnames(id_filled)[-(1:7)]
no_empty_rows<- id_filled[rowSums(is.na(id_filled[vars])) != length(vars), ]
# NA values in variables:
mean(is.na(no_empty_rows %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8658041
colnames_cleaned <- no_empty_rows %>% clean_names()

NA Values:

# NOT OK
## ten kod daje oczekiwany wynik, ale z pewnoscia da sie (trzeba) go przerobic
### czy nazwa R wziela sie od slowa "rozpacz"? 

clean_NA_in_column<-function(column){
  not_NA_count<-sum(!is.na(column))
  if (not_NA_count>=2){
    column <- na_interpolation(column, option = "linear")
    column
  }

  else if (not_NA_count==1){
    val <- first(na.omit(column))
    column[is.na(column)] <- val
    column
  }
  column
}

clean_NA <- function(input){
  res <- sapply(input, function(x) clean_NA_in_column(x))
  res
}

vars <- colnames(colnames_cleaned)[-(1:7)]
colnames_cleaned_dim<- dim(colnames_cleaned)
cleaned <- data.frame(matrix(ncol=81,nrow=0))
for (p_id in distinct(colnames_cleaned, patient_id)$patient_id){
#for (p_id in 1:10){
  patient_data <- colnames_cleaned%>%filter(patient_id==p_id)
  row_count<- dim(patient_data)[1]
  dt=data.frame(matrix(ncol=0,nrow=row_count))
  for (var in vars){
    column <- patient_data%>%select(var)
    res<-clean_NA_in_column(column)
    dt<-cbind(dt,res)
  }
  patient_data[-(1:7)] = dt
  cleaned<-rbind(cleaned, patient_data)
}
## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(var)` instead of `var` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

2.5 Prezentacja czystych danych

“Sekcja podsumowująca rozmiar zbioru i podstawowe statystyki” Szczegółowa analiza wartości atrybutów.

DT::datatable(cleaned, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))
## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html
all_dim <-dim(cleaned)
tmp <- cleaned %>% select(patient_id, outcome, gender) %>% group_by(patient_id) %>% summarise(outcome_count = mean(outcome), gender_count =  mean(gender))
## `summarise()` ungrouping output (override with `.groups` argument)
# ile pacjentow
patients_count <- length(distinct(cleaned, patient_id)$patient_id)
# ile pomiarow
measurements_count <- all_dim[1]
# ile srednio pomiarow na pacjenta
mean_measures <- measurements_count/patients_count
# ile kobiet/mezczyzn
genders <- tmp %>% count(gender_count) #Male 224, Female 151 # WYKRES
# ile umarlo/przezylo
outcomes <- tmp %>% count(outcome_count) #201 recovered, 174 died # WYKRES
# ile kolumn
columns_count <- all_dim[2]
# ile atrybutow
vars_count <- all_dim[2]-7
# ile zostalo wartosci NA
na_left <- round(100*mean(is.na(cleaned)), 0)

titles <- c("Liczba pacjentów", "Liczba pomiarów", "Średnia liczba pomiarów", "K", "M", "Smierc", "Zycie :(", "Liczba wierszy", "Liczba zmiennych", "Brakujace wartosci [%]")
values <- c(patients_count,measurements_count,mean_measures,1, 2, 0,1, columns_count,vars_count, na_left)
info_table <- data_frame(
  titles,
  values
)
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
knitr::kable(info_table)
titles values
Liczba pacjentów 361.00000
Liczba pomiarów 5606.00000
Srednia liczba pomiarów 15.52909
K 1.00000
M 2.00000
Smierc 0.00000
Zycie :( 1.00000
Liczba wierszy 81.00000
Liczba zmiennych 74.00000
Brakujace wartosci [%] 8.00000
#wykres death_bothgenders, alive_bothgenders
invisible(remove(tmp))

plot_data <- cleaned%>%select(patient_id, gender, admission_time, discharge_time, outcome) %>%group_by(patient_id) %>% summarise_all(funs(first))
## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
admission_plot <- ggplot() + 
    coord_cartesian() +
    scale_color_hue() +
    facet_wrap(~outcome) +
    layer(data=plot_data, 
          mapping=aes(
              x=admission_time, 
              y=discharge_time, colour=factor(gender)
              ), 
          stat="identity", 
          geom="point", 
          position=
              position_jitter()
    )
ggplotly(admission_plot)
## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

2.6 Analiza wartosci atrybutow

Szczegółową analizę wartości atrybutów. TBA

2.7 Korelacja miedzy danymi

“Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji.”

# OK
#TODO wiecej wynikow
var_df<-cleaned[-(1:7)]
plot_data <- corr_cross(var_df, 
max_pvalue = 0.05,
top = 30,
rm.na=TRUE)
## Returning only the top 30. You may override with the 'top' argument
## Warning in theme_lares2(legend = "top"): Font Arial Narrow is not installed, has
## other name, or can't be found
#plot(plot_data)
ggplotly(plot_data)

2.8 Zmiana [atrybutu/ow] w czasie

Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie.

timeline_plot <- ggplot() + 
    coord_cartesian() +
    scale_color_hue() +
    layer(data=cleaned, 
          mapping=aes(
              x=re_date, 
              y=hypersensitive_cardiac_troponin_i, group=patient_id
              ), 
          stat="identity", 
          geom="point", 
          position=
              position_jitter()
    )
ggplotly(timeline_plot)
timeline_plot <- ggplot(cleaned, aes(x=re_date, y=serum_chloride, colour=factor(patient_id), group=patient_id))  + geom_line() + geom_point() + facet_wrap(~outcome)
ggplotly(timeline_plot)
#timeline_data  <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose)
timeline_data  <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose, outcome, gender) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose, outcome, gender) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose), outcome=first(outcome), gender=first(gender)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose, outcome, gender)
## `summarise()` regrouping output by 'id' (override with `.groups` argument)
timeline_plot <- ggplot(timeline_data,aes(day_no,avg_hemoglobin, shape=factor(gender)))+geom_point(color="blue")+geom_point(aes(day_no,avg_glucose, shape=factor(gender)),color="black") + facet_wrap(~outcome)
ggplotly(timeline_plot)
## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.

2.9 Klasyfikator

klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji.

2.9.1 Analizę ważności atrybutów najlepszego znalezionego modelu

2.10 Dodatkowa analiza

analiza typowa dla danych klinicznych, np.:

  • regresja logistyczna wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors)
  • regresja Coxa (ang. Cox Proportional-Hazards Model).